PROJECT SETUP

PART 1 - About the data

Sample Metadata

The data being used in this workflow is taken from “Spatial Gene Expression Changes in the Mouse Heart After Base-Targeted Irradiation”

Radiation cardiotoxicity (RC) is a clinically significant adverse effect of treatment for patients with thoracic malignancies. Clinical studies in lung cancer have indicated that heart substructures are not uniformly radiosensitive, and that dose to the heart base drives RC.

This was a pilot experiment to determine whether spatial transcriptomics could provide useful insights in this context.

The data is included in the posit cloud workshop but if you are running this workshop elsewhere you can down load it from dropbox

Experimental Design:

A subset of the data is used for this workshop and it includes:

2 samples 2 conditions - irradiated heart and control 2 adjacent sections taken from each tissue block (sample) 4 sections in total

The sample metadata is summarised below:

dataDir <- "testData/results"

samples <- list.files(dataDir)


samples_short <- samples %>% 
  gsub("results_", "", .)  %>% 
  gsub("_Results", "", .)

sample_id <- samples_short %>% 
  gsub("_.*", "", .)  



condition <- c(rep("Control", 2), rep("irradiated", 2))
tissue <- c(c(rep("tissue1", 2), rep("tissue2", 2)))

sampleInfo <- data.frame(Fname = samples, 
          id = samples_short, 
          condition, 
          tissue, 
          sample_id
          )

sampleInfo
##                                         Fname                          id
## 1    results_C01_QUB_002_B1NP_Control_Results    C01_QUB_002_B1NP_Control
## 2    results_C02_QUB_002_B1NP_Control_Results    C02_QUB_002_B1NP_Control
## 3 results_C05_QUB_003_B2BP_irradiated_Results C05_QUB_003_B2BP_irradiated
## 4 results_C06_QUB_003_B2BP_irradiated_Results C06_QUB_003_B2BP_irradiated
##    condition  tissue sample_id
## 1    Control tissue1       C01
## 2    Control tissue1       C02
## 3 irradiated tissue2       C05
## 4 irradiated tissue2       C06

Spaceranger

The data is exported from the 10X pipeline Spaceranger and it conists of the aligned and quantified reads for the transcriptomics data and the processed image data.

pathToTenXOuts <- file.path(dataDir, sampleInfo$Fname, "outs")

Import the expression data

The SpatialFeatureExperiment library contains an import function for 10X data which takes the path to 10X files and loads the data in as a SpatialFeatureExperiment object

sfe <- SpatialFeatureExperiment::read10xVisiumSFE(
                                      dirs = pathToTenXOuts,
                                      samples = sampleInfo$Fname,
                                      sample_id = sampleInfo$sample_id,
                                      type = "sparse",
                                      data = "filtered",
                                      images = "hires",
                                      zero.policy=TRUE)
saveRDS(sfe, "rObjects/sfe.rds")

The import function is a bit time consuming so save time I have pre-computed a loaded data rds file which can be imported here:

sfe <- readRDS("rObjects/sfe.rds")
sfe
## class: SpatialFeatureExperiment 
## dim: 31053 14053 
## metadata(0):
## assays(1): counts
## rownames(31053): ENSMUSG00000051951 ENSMUSG00000089699 ...
##   ENSMUSG00000096730 ENSMUSG00000095742
## rowData names(1): symbol
## colnames(14053): AAACAAGTATCTCCCA-1 AAACAGAGCGACTCCT-1 ...
##   TTGTTTGTATTACACG-1-3 TTGTTTGTGTAAATTC-1-1
## colData names(4): in_tissue array_row array_col sample_id
## reducedDimNames(0):
## mainExpName: NULL
## altExpNames(0):
## spatialCoords names(2) : pxl_col_in_fullres pxl_row_in_fullres
## imgData names(4): sample_id image_id data scaleFactor
## 
## unit: full_res_image_pixel
## Geometries:
## colGeometries: spotPoly (POLYGON) 
## 
## Graphs:
## C01: col: visium
## C02: col: visium
## C05: col: visium
## C06: col: visium

SpatialFeatureExperiment Class

Built of SpatialExperiment Class (which is built on SingleCellExperiment)

Incorporates geometries and geometry operations with the sf package

SpatialExperimentFeature Class ## SP and SF

sp and sf are to classes used in Spatial Analysis. sp is the older of the 2 and uses S4 objects, it has been around for 20 years but there is a move towards sf (simple feature) that sf uses S3 classes.

We currently use both sp and sf in our package but we are currently working towards migrating to 100% sf workflow
https://github.com/r-spatial/sf/wiki/migrating

Add spot positions to sfe

This function adds the spot positions to the as the spot coordinates which are used for plotting. For more details the functions are included R/sfe_prepareSFE.R

sfe <- Spaniel::prepareSFE(sfe, sampleInfo)
## Joining with `by = join_by(sample_id)`

PART 2 - Annotating image data and compare with clustering results

Add Image annotations

Different regions of the heart were annotated using the imaging software Amira. These could also be created using Photoshop or other imaging software. Each annotation was then exported as a separate layer. The annotation images are the same size as the high resolution images exported from

For example this image shows the left ventricle in sample “C01”:

C01

Left Ventricle sample C01

The annotation_images directory

annotationDir <- "testData/annotation_images/"
annotations <- list.files(annotationDir) 
annoInfo <- data.frame(sample_id = gsub("_.*", "", annotations) ,
  domain = gsub(".*_", "", annotations) %>% 
    gsub("\\.jpg", "", .),jpg = annotations) 



annoInfo
##    sample_id          domain                     jpg
## 1        C01          centre          C01_centre.jpg
## 2        C01     left-atrium     C01_left-atrium.jpg
## 3        C01  left-ventricle  C01_left-ventricle.jpg
## 4        C01    right-atrium    C01_right-atrium.jpg
## 5        C01 right-ventricle C01_right-ventricle.jpg
## 6        C01          valves          C01_valves.jpg
## 7        C01           XR-LV           C01_XR-LV.jpg
## 8        C01           XR-RV           C01_XR-RV.jpg
## 9        C02          centre          C02_centre.jpg
## 10       C02     left-atrium     C02_left-atrium.jpg
## 11       C02  left-ventricle  C02_left-ventricle.jpg
## 12       C02    right-atrium    C02_right-atrium.jpg
## 13       C02 right-ventricle C02_right-ventricle.jpg
## 14       C02          valves          C02_valves.jpg
## 15       C02           XR-LV           C02_XR-LV.jpg
## 16       C02           XR-RV           C02_XR-RV.jpg
## 17       C05          centre          C05_centre.jpg
## 18       C05     left-atrium     C05_left-atrium.jpg
## 19       C05  left-ventricle  C05_left-ventricle.jpg
## 20       C05 right-ventricle C05_right-ventricle.jpg
## 21       C05          valves          C05_valves.jpg
## 22       C05           XR-LV           C05_XR-LV.jpg
## 23       C05           XR-RV           C05_XR-RV.jpg
## 24       C06          centre          C06_centre.jpg
## 25       C06     left-atrium     C06_left-atrium.jpg
## 26       C06  left-ventricle  C06_left-ventricle.jpg
## 27       C06    right-atrium    C06_right-atrium.jpg
## 28       C06 right-ventricle C06_right-ventricle.jpg
## 29       C06          valves          C06_valves.jpg
## 30       C06           XR-LV           C06_XR-LV.jpg
## 31       C06           XR-RV           C06_XR-RV.jpg

Add radiation field

Radiation field = black rectangle; right atrium = yellow; left atrium = blue; great vessels = red; right ventricle = purple; left ventricle = green; atrioventricular valves = orange; central fibrous body = black; irradiated right ventricle = white; irradiated left ventricle = pink.

annoInfo$radiation_field <- ifelse(annoInfo$domain %in% c("left-ventricle", 
                                                          "right-ventricle"), 
                                   yes = FALSE, no = TRUE)
annoInfo
##    sample_id          domain                     jpg radiation_field
## 1        C01          centre          C01_centre.jpg            TRUE
## 2        C01     left-atrium     C01_left-atrium.jpg            TRUE
## 3        C01  left-ventricle  C01_left-ventricle.jpg           FALSE
## 4        C01    right-atrium    C01_right-atrium.jpg            TRUE
## 5        C01 right-ventricle C01_right-ventricle.jpg           FALSE
## 6        C01          valves          C01_valves.jpg            TRUE
## 7        C01           XR-LV           C01_XR-LV.jpg            TRUE
## 8        C01           XR-RV           C01_XR-RV.jpg            TRUE
## 9        C02          centre          C02_centre.jpg            TRUE
## 10       C02     left-atrium     C02_left-atrium.jpg            TRUE
## 11       C02  left-ventricle  C02_left-ventricle.jpg           FALSE
## 12       C02    right-atrium    C02_right-atrium.jpg            TRUE
## 13       C02 right-ventricle C02_right-ventricle.jpg           FALSE
## 14       C02          valves          C02_valves.jpg            TRUE
## 15       C02           XR-LV           C02_XR-LV.jpg            TRUE
## 16       C02           XR-RV           C02_XR-RV.jpg            TRUE
## 17       C05          centre          C05_centre.jpg            TRUE
## 18       C05     left-atrium     C05_left-atrium.jpg            TRUE
## 19       C05  left-ventricle  C05_left-ventricle.jpg           FALSE
## 20       C05 right-ventricle C05_right-ventricle.jpg           FALSE
## 21       C05          valves          C05_valves.jpg            TRUE
## 22       C05           XR-LV           C05_XR-LV.jpg            TRUE
## 23       C05           XR-RV           C05_XR-RV.jpg            TRUE
## 24       C06          centre          C06_centre.jpg            TRUE
## 25       C06     left-atrium     C06_left-atrium.jpg            TRUE
## 26       C06  left-ventricle  C06_left-ventricle.jpg           FALSE
## 27       C06    right-atrium    C06_right-atrium.jpg            TRUE
## 28       C06 right-ventricle C06_right-ventricle.jpg           FALSE
## 29       C06          valves          C06_valves.jpg            TRUE
## 30       C06           XR-LV           C06_XR-LV.jpg            TRUE
## 31       C06           XR-RV           C06_XR-RV.jpg            TRUE

The annotations can be loaded in and added to the sfe object using the spaniel function findAllDomains. This creates a separate for each histological domain in the colData of the object. The combineDomains function can then be used to combine the annotations into a single column named domain.

These functions are found in: R/sfe_findDomain.R

For each of the images the functions:

  1. Load the image
  2. Clean the image to remove noise
  3. Generate the boundary coordinates of the image
  4. Convert the domain coordinates to sp (spatial polygon)
  5. Check whether spot falls inside the domain (add to metadata in sfe)
  6. Convert domain sp to sf and store within sfe object
sfe <- Spaniel::findAllDomains(sfe, annoInfo, annotationDir)
## Warning in imager::grayscale(im): Image appears to already be in grayscale mode
## Warning in imager::grayscale(im): Image appears to already be in grayscale mode
## Warning in imager::grayscale(im): Image appears to already be in grayscale mode
## Warning in imager::grayscale(im): Image appears to already be in grayscale mode
## Warning in imager::grayscale(im): Image appears to already be in grayscale mode
## Warning in imager::grayscale(im): Image appears to already be in grayscale mode
## Warning in imager::grayscale(im): Image appears to already be in grayscale mode
## Warning in imager::grayscale(im): Image appears to already be in grayscale mode
## Warning in imager::grayscale(im): Image appears to already be in grayscale mode
## Warning in imager::grayscale(im): Image appears to already be in grayscale mode
## Warning in imager::grayscale(im): Image appears to already be in grayscale mode
## Warning in imager::grayscale(im): Image appears to already be in grayscale mode
## Warning in imager::grayscale(im): Image appears to already be in grayscale mode
## Warning in imager::grayscale(im): Image appears to already be in grayscale mode
## Warning in imager::grayscale(im): Image appears to already be in grayscale mode
## Warning in imager::grayscale(im): Image appears to already be in grayscale mode
## Warning in imager::grayscale(im): Image appears to already be in grayscale mode
## Warning in imager::grayscale(im): Image appears to already be in grayscale mode
## Warning in imager::grayscale(im): Image appears to already be in grayscale mode
## Warning in imager::grayscale(im): Image appears to already be in grayscale mode
## Warning in imager::grayscale(im): Image appears to already be in grayscale mode
## Warning in imager::grayscale(im): Image appears to already be in grayscale mode
## Warning in imager::grayscale(im): Image appears to already be in grayscale mode
## Warning in imager::grayscale(im): Image appears to already be in grayscale mode
## Warning in imager::grayscale(im): Image appears to already be in grayscale mode
## Warning in imager::grayscale(im): Image appears to already be in grayscale mode
## Warning in imager::grayscale(im): Image appears to already be in grayscale mode
## Warning in imager::grayscale(im): Image appears to already be in grayscale mode
## Warning in imager::grayscale(im): Image appears to already be in grayscale mode
## Warning in imager::grayscale(im): Image appears to already be in grayscale mode
## Warning in imager::grayscale(im): Image appears to already be in grayscale mode
## create a column containing all domains
sfe <- combineDomains(sfe, unique(annoInfo$domain))

sfe$domain

## add radiation field
sfe$radiation_field <- colData(sfe) %>% 
  data.frame() %>% 
  left_join(annoInfo) %>% 
  pull(radiation_field)
## Joining with `by = join_by(sample_id, domain)`

The domains can now be visualised on the histological data

p1 <-  Spaniel::spanielPlot_SFE(sfe,  
                 sample_id = "C01", 
                 plot_feat = "domain", ptSizeMin = 0, 
                 ptSizeMax = 2) 
## [1] "got coordinates"
## [1] "metadata"
## [1] "added plot feature"
## [1] "discrete"
## SpatRaster resampled to ncells = 501264
p2 <-  Spaniel::spanielPlot_SFE(sfe,  
                 sample_id = "C01", 
                 plot_feat = "radiation_field", ptSizeMin = 0, 
                 ptSizeMax = 2) 
## [1] "got coordinates"
## [1] "metadata"
## [1] "added plot feature"
## [1] "discrete"
## SpatRaster resampled to ncells = 501264
p1 + p2

Sample QC

Now that the image annotations are fully loaded in the transcriptomics data can be analysed.

QC metrics about each spot can be calculated using the addPerCellQCMetrics function from scuttle package. By default the total genes per spot are calculated. Additional metrics such as the percentage of mitochondrial genes can be calculated. As this data is from mouse mitochondrial genes start with the prefix “mt-” and can be extracted based on pattern matching.

is_mt <- str_detect(rowData(sfe)$symbol, "^mt-")
sfe <- scuttle::addPerCellQCMetrics(sfe, subsets = list(mito = is_mt))

Spots where no genes are detected can be removed from the remainder of the analysis.

filter <- sfe$detected > 0 & sfe$in_tissue
sfe <- sfe[, filter]

Split the data into control and irradiated

Although we could now progress with the analysis using all samples it is often advisable to analyse each tissue sample or condition separately to begin with. This allows you to identify any section or sample which might be an outlier and should be removed from downstream analysis. For this workshop we example we will split the samples by condition:

sfe_control <- sfe[,sfe$condition == "Control"]
sfe_irradiated <- sfe[,sfe$condition == "irradiated"]

Cluster analysis

The filtered data can then be then be normalised and dimension reduction steps applied and the spots can be clustered according to transcriptomic similarity. This analysis uses methods taken from the bioconductor Orchestrating Single-Cell Analysis https://bioconductor.org/books/release/OSCA/ work flows. Wrapper function (R/sfe_helper_functions.R) to perform this analysis are implemented in Spaniel for convenience. However for a more comprehensive clustering analysis we recommend following the reading the Scran and Bluster documentation which describes each individual step in the analysis in more detail.

In this example we use a Leidan clustering algorithm as implemented in scran. Leidan is a community detection method. It is important to note that the K parameter is the number of nearest neighbors to use when creating the k nearest neighbor graph. k is related to the resolution of the clustering result, a bigger k will result in lower resolution and vice versa.

sfe_control <- sfe_control %>% 
  Spaniel::preProcess() %>%
  Spaniel::clustCells(K = 100, cluster_function = "leiden")
sfe_irradiated <- sfe_irradiated %>% 
   Spaniel::preProcess() %>% 
   Spaniel::clustCells(K = 100, cluster_function = "leiden")

Control data

Now we can look at the results of the data both by visualising the clusters on a UMAP plot using plot functions from the scater package and within the spatial context of the tissue section using the spanielPlot function (R/sfe_PlotFunctions.R).

p1_control <- Spaniel::spanielPlot_SFE(sfe_control,  
                 sample_id = "C01", 
                 plot_feat = "domain", ptSizeMin = 0, 
                 ptSizeMax = 2) 
## [1] "got coordinates"
## [1] "metadata"
## [1] "added plot feature"
## [1] "discrete"
## SpatRaster resampled to ncells = 500472
p2_control <- Spaniel::spanielPlot_SFE(sfe_control,  
                 sample_id = "C01", 
                 plot_feat = "clust_leiden_K_100", 
                 ptSizeMin = 0, 
                 ptSizeMax = 2) 
## [1] "got coordinates"
## [1] "metadata"
## [1] "added plot feature"
## [1] "discrete"
## SpatRaster resampled to ncells = 500472
p1_control + p2_control

irradiated data

scater::plotUMAP(sfe_irradiated, colour_by = "domain")

p1_irradiated <- Spaniel::spanielPlot_SFE(sfe_irradiated,  
                 sample_id = "C06", 
                 plot_feat = "domain", ptSizeMin = 0, 
                 ptSizeMax = 2) 
## [1] "got coordinates"
## [1] "metadata"
## [1] "added plot feature"
## [1] "discrete"
## SpatRaster resampled to ncells = 501028
p2_irradiated <- Spaniel::spanielPlot_SFE(sfe_irradiated,  
                 sample_id = "C06", 
                 plot_feat = "clust_leiden_K_100", 
                 ptSizeMin = 0, 
                 ptSizeMax = 2) 
## [1] "got coordinates"
## [1] "metadata"
## [1] "added plot feature"
## [1] "discrete"
## SpatRaster resampled to ncells = 501028
p1_irradiated + p2_irradiated

Trying a lower resolution

sfe_control <- sfe_control %>% 
  Spaniel::clustCells(K = 75, cluster_function = "leiden")
sfe_irradiated <- sfe_irradiated %>% 
   Spaniel::clustCells(K = 75, cluster_function = "leiden")
p3_control <- Spaniel::spanielPlot_SFE(sfe_control,  
                 sample_id = "C01", 
                 plot_feat = "clust_leiden_K_75", 
                 ptSizeMin = 0, 
                 ptSizeMax = 2) 
## [1] "got coordinates"
## [1] "metadata"
## [1] "added plot feature"
## [1] "discrete"
## SpatRaster resampled to ncells = 500472
p3_irradiated <- Spaniel::spanielPlot_SFE(sfe_irradiated,  
                 sample_id = "C06", 
                 plot_feat = "clust_leiden_K_75", 
                 ptSizeMin = 0, 
                 ptSizeMax = 2)
## [1] "got coordinates"
## [1] "metadata"
## [1] "added plot feature"
## [1] "discrete"
## SpatRaster resampled to ncells = 501028
p1_control + p2_control + p3_control

p1_irradiated + p2_irradiated + p3_irradiated

Save the processed data

saveRDS(sfe_control, "rObjects/wrapped_sfe_control.rds")
saveRDS(sfe_irradiated, "rObjects/wrapped_sfe_irradiated.rds")

PART 3 - Share analysed project

Prepare the data to share

Finally we prepare the data so that it can be loaded into the Shiny Spaniel app. To demonstrate this in the interest of time are going to look at a subset of the data.

sfe_single <- sfe_control[,sfe_control$sample_id == "C01"]

The app allows us to look at QC features, and browse between different clustering solutions as well as different sample conditions. Before we share our analysis with the rest of the team we can mark these columns in the metadata using the spaniel helper function “labelSFE”.

We can check what we have in the metadata using the following function:

colnames(colData(sfe_single))
##  [1] "in_tissue"             "array_row"             "array_col"            
##  [4] "sample_id"             "Fname"                 "id"                   
##  [7] "condition"             "tissue"                "centre"               
## [10] "left-atrium"           "left-ventricle"        "right-atrium"         
## [13] "right-ventricle"       "valves"                "XR-LV"                
## [16] "XR-RV"                 "domain"                "radiation_field"      
## [19] "sum"                   "detected"              "subsets_mito_sum"     
## [22] "subsets_mito_detected" "subsets_mito_percent"  "total"                
## [25] "sizeFactor"            "clust_leiden_K_100"    "clust_leiden_K_75"

Then we can define the columns that we want to label in our data for QC plots, groups (clusters or anatomical domains, histological regions) and the experimental conditions in our experiment.

grps <- c("clust_leiden_K_100", "clust_leiden_K_75", "domain", "radiation_field")
qc <- c("total", "sum", "subsets_mito_percent", "in_tissue")
exp_cond <- c("condition")

sfe_single <- Spaniel::labelSFE(sfe_single, qc_cols = qc, grp_cols = grps, exp_cond_cols = exp_cond)

As you can see the colData columns are now updated:

colnames(colData(sfe_single))
##  [1] "qc_in_tissue"            "array_row"              
##  [3] "array_col"               "sample_id"              
##  [5] "Fname"                   "id"                     
##  [7] "expCond_condition"       "tissue"                 
##  [9] "centre"                  "left-atrium"            
## [11] "left-ventricle"          "right-atrium"           
## [13] "right-ventricle"         "valves"                 
## [15] "XR-LV"                   "XR-RV"                  
## [17] "grp_domain"              "grp_radiation_field"    
## [19] "qc_sum"                  "detected"               
## [21] "subsets_mito_sum"        "subsets_mito_detected"  
## [23] "qc_subsets_mito_percent" "qc_total"               
## [25] "sizeFactor"              "grp_clust_leiden_K_100" 
## [27] "grp_clust_leiden_K_75"
saveRDS(sfe_single, "rObjects/sfe_single.rds")

Run Shiny Spaniel!

Now we can launch the shiny app directly from rstudio

Run shiny spaniel

Spaniel::runShinySpaniel_SFE()